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ABSTRACT 

It is now established that several of the Local Group dwarf Spheroidal galaxies (dSphs) 
have large mass-to-light ratios. We consider the possibility that the dark matter in the 
haloes of dSphs is composed of massive black holes with masses in the range 10 5 M 
to 10 7 M Q . We use direct iV-body simulations to determine the long term evolution 
of a two-component dSph composed of a pressure-supported stellar population within 
a black hole dominated halo. The black holes are initially distributed according to a 
Navarro, Frenk & White profile. For black hole masses between 1O 5 M0 and 1O 6 M0, 
the dark matter halo evolves towards a shallower inner profile in less than a Hubble 
time. This suggests that black holes in this mass range might provide an explanation 
for the origin of the dark matter cores inferred from observations of Low Surface 
Brightness galaxy rotation curves. We compare the simulated evolution of the stellar 
population with observed data for the Draco dSph and find that dynamical heating 
generally leads to the rapid dispersal of the stellar population to large radii. The 
dependence of the heating rate on the black hole mass is determined, and an upper 
limit of 10 5 Mq is placed on the individual black holes comprising the dark matter 
halo of Draco, if its present stellar distribution is representative of the initial one. 
We also present a simple scaling argument which demonstrates that the dynamical 
heating of an initially compact, though not self-gravitating, stellar distribution might 
produce a remnant distribution similar in extent to Draco after 10 Gyr, even for black 
hole masses somewhat in excess of 10 5 Mq. 

Key words: dark matter — galaxies: dwarf, haloes — methods: TV-body simulations 



1 INTRODUCTION 

The existence of dark matter in the Universe, which has 
an abundance six times that of baryonic matter, is inferred 
from sources as diverse as the disk rotation curves of spi- 
ral galaxies, gravitational lensing studies of galaxy clusters, 
the very high mass-to-light ratios of nearby dwarf galaxies 
and analysis of Cosmic Microwave Background fluctuations. 
However, the actual composit ion of the dark matter in th e 
Universe is as yet unknown JOstriker fc Steinhardtl feoOS). 
One method by which dark matter may be detected and its 
properties probed is by means of substructu r e lensing (e.g . 
iMetcalf fc Madaull200ll : iLi fc Ostrikerll2002l lKeetonll2003) . 
Dark matter studies in the Local Group can also yield con- 
straints on dark matter models via the m ass density as in- 
ferred from kinematic measurements (e.g. IWilkinson et alJ 
120021) . 

The standard cold dark matter (CDM) cosmology, 
though popular and well-tested on large scales, has some un- 
resolved issues. The two best known are the over-abundance 
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of low-mass satellite haloes identified at the present epoch 
in cosmological simulations relative to the number of ob - 
served Local Group satellites (e.g. iMoore et alJ Il999l) 
and the much-debated inner profile of CD M haloes (e.g. 
Ide Blok fc Bosmal2002llSwaters et al.l2005t) . The latter may 
be due to resolution issues in the simulations, where a 
range of inner densi ty slopes are found by various au- 
thors (see ISutoll2002l for a concise summary). Low mass 
haloes may also have intri nsically softer inner profiles than 
more massive systems jRicottil l200otl . Observationally, 
there is considerable evidence that massive galaxies also 
do not have cusped haloes ; but instead contain roughly 
uniform density cores _(e^g. iBinnev fc Evansll200lt Salucci 
l200ltlBorriello et all20ol) . The issue of the missing satellite 
haloes may be due to the fact that some of the satellite galax- 
ies are very diffuse; a few bound stellar systems may there- 
fore be undet ected simply due to their low surface brightness 
iMateolll998l) . It has also been suggested that the satellite 
crisis can be resolved if all the Local Group dwarf spheroidal 
galaxies (dSphs) lie in significantly more massive haloes than 
those implied by previous estimates, whi ch were based only 
on their central velocity dispersions (e.g. IStoehr et alJl2003 : 
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lHavashi et al However, differences between the radial 

distribution of massive satellite haloes about their parent 
galaxies in simulations and the observed distribution of the 
Milky Way dSphs suggest t hat this cannot be the complete 
solution to this problem (e.g. lDe Lucia et al.l20o4) . Another 
possibility is that there may exist dark matter haloes which 
do not harbour a stellar population. 

In order to address these issues, some alternatives to 
CDM have been proposed: these incl ude warm dark matter 
and self-interacting dark matte r (e.g. ISpergel fc Steinhardd 
1200(1 iGnedin fc Ostrikerll200 lfl . An interesting complete al- 
ternative to dark matter has also been suggested for the 
observed dwarf Spheroidals (dSphs) whereby special phase- 
space characteristics permit the long-lived remnant of a 
tidally disrupted sa tellite to assume the ap pearance of a 
dSph in projection Jki CSSCll fc Kroupalll998l) ; however, dif- 
ficulties in reproducin g the data on Draco have recently been 
noted in this model (Klesscn et al. 20q3). 

There has been considerable interest in cosmological 
CDM simulations in the past decade, with a variety of halo 
profile fits being presented by different authors. Following a 
set o f large cosmological simulations, Na varro et al] (1996, 
Il997f) proposed a universal density profile to accommodate 
virialised CDM haloes of all masses, which has become com- 
monly known as the NFW profile. The density in the in- 
ner region is cusped as r , whereas at larger radii it is 
proportional to r -3 , with the location of the transition be- 
ing defined by a scale radius used to characterise the pro- 
file. Several authors have subsequently confirmed the fit 
provided by the NFW profile, although others have also 
presented steeper inner slopes for the halo density (e.g. 
[ Fukushige fc Makinoll200ll : iGhigna et alll200(l IColin et al] 
2003). A generalised density law has therefore been sug- 
gested, for which the original NFW profile is a special case, 
in order to account for such findings (e.g. Zhac3 ll996r) . As 
long as the dark matter is cold and collisionless, the halo 
profile would not be expected to depend on the nature of the 
dark matter and, more specifically, it should be independent 
of the dark matter particle mass. This remains true unless 
we consider particles whose masses are much greater than 
a solar mass, for which two-body dynamical effects would 
become relevant. Such phenomena would first be observed 
in the lowest mass, lowest velocity dispersion systems such 
as the dSphs. 

Much interest in recent years has been focused on the 
study of the Local Group satellites, not least because their 
relative proximity h as led to the increased availability of 
high quality data lllrwin fc Hatzidimitrioul Il995l iMateol 
1998), which place stricter constraints on the nature of 
the dark matter that appears to dominate their gravitat- 
ing mass. These systems have also received m uch attention 
on the modelling and simu lation front (e.g. IStoehr et al] 
120021: 1 Wilkinson et al]l2002h . Amongst those with the high- 
est mass-to-light ratios is the Draco dSph galaxy. At a he- 
liocentric distance of approximately 80 kpc, Draco has been 
studied extensively and a wealth of kinematic data exists; 
in particular, measurem ents of stellar velocities out to large 
radii are now available llKlevna et all2002t IWilkinson et al] 
12004) . Given that within the region probed by the stellar 
distribution it also appears to have suff ered relatively little 
from the tidal effects of the Milky Way iOd enkirchc n et al] 
l2001al : iKlessen et all 120031: IWilkinson et al]|2004 . unlike 



some of the other Local Group satellites such as the Sagit- 
tarius dwarf (e.g. iHelmi fc WhitdEoOll) . Draco is an ideal 
candidate with which to study the possible properties of 
dark matter. Draco is also known to contain an old stellar 
population with an age of more than 10 Gyr; any model 
of this galaxy must therefore remain stable for at least this 
length of time in order to be compatible with this observa- 
tion. 

One conceivable candidate for the dark matter, either 
as the sole or major component, is an ensemble of massive 
black holes. These could be primeval, or may have grown 
by accretion and mergers from a prime val population, or 
have origins in later non- linear structures. lLacev fc Ostrikerl 
(1985: hereafter L085) investigated the way in which a halo 
of massive black holes would affect the stellar dynamics in 
the Milky Way and obtained results that are generally con- 
sistent with observations. Our present study investigates the 
dynamical effects of a dark matter halo composed purely of 
primordial massive black holes on the evolution of the stel- 
lar component of a dwarf galaxy. It was shown in L085 that 
for a dSph such as Draco, an upper limit for the individual 
black hole mass would be 2 x 10 6 Mq, based upon stabil- 
ity arguments for both the stellar population an d the dark 
matte r halo forming the galaxy. More recently, iMao et ail 
(2004) have suggested that massive black holes with masses 
in the range 10 5 — 10 6 Mq may offer a viable explanation for 
the observed image flux ratios in strong gravitational lens 
systems. We use direct iV-body simulations to trace the evo- 
lution of a model dSph whose initial structural parameters 
are chosen such that they are compatible with present-day 
observations of Draco. The stellar population is treated as 
a cluster of tracer particles, whose evolution is dependent 
on the potential of the dark halo. The distribution of the 
halo particles is chosen to follow an NFW profile, whilst the 
stellar population is given an initial distribution which is 
described by a Plummer law. We assume the dark matter of 
the halo to consist entirely of massive black holes, formed 
prior to the formation of the stellar population. Their masses 
are chosen to cover a range whose maximum value has, for 
generality, been taken to be larger than the constraints im- 
posed by L085; in each simulation, all black holes have the 
same mass. 

In Section |21 we discuss the generation of our simulated 
models. We present the results in Section [3] along with a 
discussion in Section |3] In Section |S] we discuss alternative 
initial conditions for the stellar distribution of dSphs and 
their implications for black hole haloes. Finally we present 
a summary and our main conclusions in Section HJ 



2 SIMULATIONS 

A summary of the simulations performed is given in Ta- 
ble The simulation s were carried out using the NBODY2 
code dAarsethll200lL modified to incorporate tracer par- 
ticles orbiting within the underlying potential of the dark 
matter halo to represent the observed stars. The mass of 
the system was assumed to be derived entirely from the 
dark matter. The black hole mass was varied in different 
runs over the range 10 5 -10 7 M a . We adopt the standard - 
ised iV-body units described by iHeggie fc Mathieul (.1986): 
G — 1, Mtot = 1, where Mtot is the total self-gravitating 
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Table 1. Model parameters: Column 1 gives the simulation number; columns 2-4 give the total number of black holes, the mass of a 
single black hole and the total halo mass respectively; column 5 gives the total number of tracers; column 6 gives an estimate of the 
initial half-mass relaxation time t T ^ given by equation 1111 : column 7 gives the corresponding relaxation time for the central 5 percent 
(by mass) of the black hole distribution; columns 8-10 give the halo scale length, density scale and concentration respectively; column 
11 gives the length conversion from JV-body units to physical distances in kpc; column 12 gives the gravitational softening e (in iV-body 
units) used. 
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mass, and the initial total energy Eq = —1/4. The inclu- 
sion of a gravitational softening parameter e, or softening 
length, ensures that force singularities do not occur as the 
separation of particles tends to zero. A value of e = 0.02 
in iV-body units, chosen for computational expediency, was 
used for the majority of our simulations, corresponding to 
a physical length of 86 pc < e < 144 pc depending on the 
model. Model 5 was performed as a repeat of model 6, but 
with the softening length reduced by a factor of 100, in or- 
der to verify that we had not missed any significant physical 
effects as a result of using relatively large softening in the 
other models. The number of tracers at the start of the simu- 
lation was 10,000 for models 2-4, 8 and 11, and 100,000 for 
the remaining models. All simulations were performed on 
the Sun Workstation Cluster at the Institute of Astronomy 
and allowed to run to at least 20 Gyr. 



2.1 Initial Conditions 

The NFW density profile jNavarro et al.lll99fj fl997h . pro- 
posed as a universal fitting formula for CDM haloes which 
form in a hierarchical clustering scenario, has the form: 



p(r) = 



PsT s 



r(r + r 3 ) 2 



(1) 



where the parameters p s and r s used to characterise the 
NFW profile are the density scale and scale radius respec- 
tively. 

The black holes in our simulations are initially dis- 
tributed according to the NFW profile. The allocation of the 
initial particle velocities is performed by numerically eval- 
uating the velocity dispersion using th e Jeans equation for 
a sph erically symmetric system (e.g. iBinnev fc Tremaind 
1987), 



1 d(va 
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d$ 
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where v(r) is the halo density profile, cr(r) is the one dimen- 
sional velocity dispersion for the halo and the underlying 
halo potential is $(r). We have assumed velocity isotropy 
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Figure 1. The velocity dispersion a (solid curve) and circular 
velocity v c (dotted curve) of the halo particles as a function of 
radius. The initial velocities of the massive particles are generated 
from this distribution for all models except models 7 and 10 in 
which the halo scale length r B is larger. 



for the halo particles, which gives a satisfactory fit to the 
results observed in cosmological simulations. In this case, 
cr(r) has the form 

a 2 = 4nGp B r 2 (s max - s x )x (1 + x) 2 , (3) 
where x is the scaled radius r/r s and s x is given by 
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dx' . 



(4) 



The value of s ma x was evaluated by numerically integrating 
up to a sufficiently large radius. The form of the one dimen- 
sional dispersion so derived is shown in Fig. which also 
shows the circular velocity profile for the NFW halo. 

For each particle, we randomly generate a velocity from 
a Maxwellian velocity distribution with velocity dispersion 
calculated using equation ^ at the position of the particle. 
While at each place the total energy (gravitational and ki- 
netic) of the local particles is correct and the Jeans equation 
is sati sfied, we note that recent work bv iKazantzidis et alJ 
(2004) has found that the assumption of a Maxwellian ve- 
locity distribution can lead to a re-adjustment and short- 
term spurious evolution in iV-body simulations, due to the 
fact that the true velocity distribution is not an exact 
Maxwellian. However, this evolution does not have a sig- 
nificant impact on the simulations we present in this paper. 
The initial re-adjustment takes place on a short time scale 
(the halo crossing time) while the effects we are studying 
take place on the much longer halo relaxation time scale. 

We note in passing that it is possible to relate the max- 
ima of the rotation curve u ma x and velocity dispersion a max 
of an NFW halo, and the corresponding radii ri and T2, to 
p s and r s via 



«L X = 2.7i7G Ps r s 2 
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1.185G p s rg 



r 2 = 0.763r s 



(•>) 
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The numerical constants were obtained by direct evalua- 
tion of the circular speed and velocity dispersion profiles. 
As Fig. also shows, the maxima of the two profiles do not 
occur at the same radius. 

We represent the dSph stellar population by tracer par- 
ticles. This is justified on the basis that several of the dwarf 
satellites of the Milky Way are observed to have very high 
mass-to-light ratios; in particular, Draco is known to have 
an ext remely large mass-t o-light ratio of 350-1000 in solar 
units (iKlevna et aljl200lT) . These systems therefore appear 
to be very dark matter dominated, with their stellar popu- 
lations residing in massive dark haloes. In our simulations, 
the tracers are initially distributed according to a Plummer 
profile, which is known to provide a good fit to the present 
light distribution in the inner regions of certain dSphs, in- 
cluding Draco. The form of the three dimensional Pl ummerl 
luminosity distribution v(r) is given by 



VQTp 



+ r- 
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(7) 



where the scale length rp corresponds to the projected half- 
light radius of the model. For Draco, rp is approximately 
232 pc, which implies a core radius (the radius at which the 
surface brightness falls to half its central value) of 150 pc. 
The tracer velocities are again obtained by means of the 
Jeans equation ©. As before, we assume velocity isotropy, 
as suggested by th e observed stellar kinematics of Draco 
jKlevnaet al.ll2002l) . 

The halo parameters for the majority of our simulations 
are chose n to reflect the obse rvational estimates of the mass 
of Draco iKlevna et al]l200ll) . The halo mass is constrained 



by requiring that the value of M 2 oo, the mass contained 
within a radius r2oo which encompasses a mean overdensity 
of 200 times the critical density, yields a mass within 3rp 
which is consistent with that of Draco. The scale length of 
the halo is given by r s — T2oo /c, where the 'co ncentration' c 
of the halo is set using the scaling relation of iBullock et alJ 
feOOll) : 
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where M* = 2.14 x 10 M© is the typical collapsing mass 
at the present epoch. The value of r2oo is then given by 

1/3 

= I — i , (9) 
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where p cr it = 3H 2 /8itG is the critical density for closure of 
the Universe; we take H = 70kms _1 Mpc -1 . The relation 
between the concentration and the other parameters is 
made complete by the alternative description of the density 
scale as p s = S c p cr it where 
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200 



3 [ln(l + . 



c/(l 



(10) 



The specification of the halo mass M200 thus allows all other 
parameters to be deduced using the above relations. While 
most of our simulated haloes have similar masses to that of 
Draco, we also include a number of haloes with larger scale 
radii. This extends our modelling to include Low Surface 
Brightness (LSB) galaxies. 

2.2 Relaxation Time Scale 

We expect the time scale on which the halo profile evolves 
to be related to t r h, the relaxation time at the half-mass 
radius of the initial black hole distribution. The relaxation 
time scale is essentially the time scale over which particles 
lose memory of their initial energies due to mutual encoun- 
ters. The half-mass relaxation time for each model, listed 
in Table U is estimated using the definition of t T h as given 
by Bi nnev fc Tremaind 1119871 chapter 8) : 



trh — 



2.1 x 10V ( 10 6 M 
ln(0.4A) 



M 



10 9 M 



1 kpc 



fll) 



where M and N are the total mass and number of black 
holes in the system respectively, m is the mass of an individ- 
ual black hole and Vh is the half-mass radius. We note that 
the derivation of equation 11 H assumes a Plummer model 
and is not precisely valid for an NFW profile for which the 
ratio of the half-mass radius to the gravitational radius is 
about 15 percent larger than in the Plummer case. Thus, 
this corresponds to a deviation in i r h of less than 5% for our 
models. However, equation 111! provides a good indication 
of the typical relaxation time and we have chosen therefore 
to keep it in its more familiar form. 

The dependence of relaxation time on black hole mass 
demonstrates the conservative nature of our choice of ini- 
tial parameters for the simulated haloes. Given that mas- 
sive black holes are, by definition, manifestations of Cold 
Dark Matter, one might expect their generic density distri- 
bution to be similar to an NFW profile. However, this is 
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Figure 2. The structural evolution of the halo in model 1 as 
depicted by the Lagrange radii in Af-body units. From the low- 
est curve upwards, the radii shown correspond to the following 
halo mass percentiles: 7.5, 10, 15, 20, 50, 75 and 85 percent. The 
cluster is seen to retain its initial structure over the course of the 
simulation. 



not necessarily the case. For example, if the dark matter 
consists entirely of 10 5 Mq black holes, then a primordial 
halo of total mass 10 7 Mq contains only 100 particles. As 
we will demonstrate in the next section, a black hole halo 
with an NFW density profile develops a central core via two- 
body processes on a time scale of i r h . lOstriker fc Gnedinl 
( 1996) estimate that 10 7 Mq haloes become non-linear at 
a redshift of about z = 13 (see their Figure 1), at which 
time their virial radius can be shown to be about 0.3 kpc. 
According to equation 1111 . the relaxation time for such a 
halo t t ,7 is then about 100 Myr. If this time scale is much 
less than the difference between the times at which haloes 
of 10 7 M and 10 9 M become non-linear, then we would 
expect 10 9 Mq haloes to be formed from cored sub-haloes. 
In this case the resulting 10 9 Mq haloes would have density 
profiles which are significantly less cusped than NFW, which 
would agree with observations of LSB galaxies. In fact, the 
10 9 Mq haloes become non-linear approximately t r ,7 later 
than 10 7 M haloes and so it is less clear what the form of 
the resulting density profile should be. However, the most 
unfavourable starting point is the NFW profile and our sim- 
ulations therefore constitute the strictest test of the ability 
of black hole dark matter to produce cored galaxy haloes. 



Table 2. Halo parameters at T = and T = 15 Gyr: Column 1 
gives the simulation number; column 2 gives the logarithmic slope 
of the inner halo at T = for those models in which there are 
sufficient numbers of black holes to give a statistically meaningful 
estimate; column 3 gives the radius r max at which the initial 
circular velocity profile achieves its maximum value; column 4 
gives the maximum value of the initial circular velocity 
columns 5—7: as for columns 2-4 but for the haloes at T = 15 
Gyr. 
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3 RESULTS 

3.1 Halo Evolution 

The stability of a cluster of massive black holes with an 
NFW profile (and no tracers) was initially assessed in or- 
der to determine the validity of our assumption that such a 
cluster could act as a halo enveloping a stellar population 
(model 1). As expected from the estimate of in Table Q 
we find that the cluster retains its initial structure over the 
duration of the simulation (Fig. |5J. The inner regions are 
not reproduced in this plot: their evolution will be discussed 
in detail below. 

The evolution of the circular velocity profiles for mod- 
els 2-4, 7, 10 and 11 is shown in Fig.|3] Tabled summarises 
the properties of these curves. As will be discussed in Sec- 
tion t ne initial evolution leads to the destruction of the 
NFW cusp and the formation of a shallower core. This is 
followed by slow evolution towards core collapse - the time 
scale for this evolution varies according to the mass of the 
black holes. The inner slope of the velocity profile becomes 
steeper for models 7, 10 and 11, indicating a transition from 
the initial cusp of the NFW profile to a more uniform den- 
sity core in the central region of the halo. The evolution of 
the mass within 500 pc of the halo centre as a function of the 
initial half-mass relaxation time is shown in Fig.|l]for mod- 
els 2-4 and 11. Models 2 and 3 display a steady increase in 
the mean density within 500 pc while in model 4, which has 
better time resolution (at early times) due to its longer t T h, 
the density first decreases on a time scale of t r h- This initial 
decrease in the mass within 500 pc is followed by a slow in- 
crease as the system evolves towards core collapse. In model 
11, for which 20 Gyr corresponds to less than one relaxation 
time, the c entral density is still decreasing at the end of the 
simulation. lHavashi et all J2003T) observe similar evolution in 
their simulations of CDM haloes. They find that this varia- 
tion of the central density profile occurs on the local escape 
time scale which is approximately 136 times the local (i.e. 
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Figure 4. The evolution of the number of black holes within 
500 pc of the halo centre for models 2, 3, 4 and 11. The time 
axes of the plots have been scaled individually according to the 
half-mass relaxation time for each model. 



Figure 3. Halo rotation curves for models 2-4, 7, 10 and 11 
shown at times T = (black), 10 (green) and 20 (red) Gyr. The 
plot for model 4 also shows the circular velocities at 1 (blue) and 
2 (magenta) Gyr. In the panel for model 7, the labelled lines indi- 
cate the expected slope for an NFW and a cored density profile. 



central) relaxation time scale. As we will discuss later, we 
find that the evolution proceeds with a characteristic time 
scale comparable to , which is sim ilar in magnitude to t he 
evolutionary time scale observed hv lHavashi et all lj200.^ . 

Fig. [3] shows that, after 20 Gyr, the evolution towards 
core collapse is well advanced only for models 2 and 3. At 
the opposite end of the mass scale, model 11 has only just 
developed a core by about 10 Gyr. For model 4, we note 
that although the inner profile begins to steepen again after 
10 Gyr of evolution, the circular speed curve after 20 Gyr 
indicates that the inner regions remain less cusped than the 
initial conditions. We conclude that core collapse in black 
hole haloes will not occur within a Hubble time unless the 
black hole mass is greater than 10 6 ' 5 Mq. Most of the models 
with black hole masses below this value have the attractive 
feature that the slope of the cusp decreases with time, with 
the final approximate inner slope being shown in Table|5|for 
those models for which a statistically meaningful estimate 
can be made. Thus, NFW haloes composed of black holes 
with masses in the range 10 s Mq to 10 6 ' 5 Mq can develop 
central cores in less than a Hubble time. This may provide 
an explanation for the origin of the apparent cores in the 
density profiles of LSB galaxies. 



3.2 Evolution of the Stellar Population 

The evolution of the tracer density profile for some of our 
models is shown in Fig. where we have assumed each star 
to have a physical mass of 1 Mq . In evaluating the projected 
densities, we have taken radial bins containing 100 tracer 
particles in each of the three projected (Cartesian) planes. 
An average value was then calculated for each bin over all 
three planes. We adopt a similar approach in evaluating the 
line-of-sight velocity dispersions shown in Fig. [(J The tracer 
profile becomes extended over the course of the simulation as 
the particles acquire energy - the heating process is apparent 
from Fig. |S| where the dispersion is clearly greater for later 
times. From the evolution of the density in projection, it is 
straightforward to derive the behaviour of the core radius 
as a function of time; this is shown in Fig. [7| for models 
4 and 6-10. As Fig. |7| shows, in all cases the core radius 
rapidly increases to values which are inconsistent with the 
observed value of 150 pc for Draco. Complementary to this 
is Fig. |H] which shows the decline of the normalized tracer 
count within 500 pc as a function of time for several models. 
Unsurprisingly, reduced softening is seen to increase the rate 
at which the stars are heated; the same effect is seen if we 
increase the black hole mass. For models utilising 10 6 Mq 
black holes, an increase in the halo scale length reduces the 
rate of heating, whereas for smaller mass black holes, the 
effect is negligible. As expected, simply increasing the size of 
the stellar population by a factor of 10 produces no effect as 
regards the tracer evolution, except that we obtain a better 
statistical description of the profile. In Section 14.41 we will 
discuss how a comparison of the core radii of the simulated 
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Figure 5. Projected tracer density plots for models 2, 3, 4 and 8 
shown at times T = (open circles), 10 (filled triangles) and 20 
(open squares) Gyr. The plot for model 4 also shows the densities 
at sampling times of 1 and 2 Gyr with filled circles and open 
triangles respectively. 



Figure 7. The tracer core radius shown as function of time. The 
upper curves show models 4 (solid), 6 (dotted) and 7 (dashed); the 
lower curves refer to models 8 (solid), 9 (dotted) and 10 (dashed). 
Note that the observed core radius of Draco is ~150pc. 
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Figure 6. Line-of-sight velocity dispersion plots for models 2, 3, 
4 and 8 shown at times T = (black), 10 (green) and 20 (red) 
Gyr. 



tracer populations with that of Draco can be used to place a 
constraint on the mass of individual black holes in the halo 
of Draco. 

The majority of our simulated systems produced no 
tracer-black hole binaries, although we observe steady black 
hole-black hole binary formation at all times for all of our 
models. Tracer-black hole binaries were, however, observed 
in models 6, 7, 9 and 10, where the first two models produced 
two binaries each during the course of the simulation, whilst 
in the latter two models, six and thirteen tracer-black hole 
binaries were observed respectively. Of these, both models 
showed one long-lived binary each, corresponding to a bound 
system lasting up to 340 and 330 Myr respectively. Binaries 
composed of two black holes were observed frequently and 
consistently in all simulations, with typically 2-4 transient 
binaries per sampling time, which corresponds to a stable 
bound system lasting for approximately 170 Myr (depend- 
ing on the model) ; occasionally a binary remained intact for 
several Gyr. Simply reducing the softening (but not down 
to zero) did not appear to increase the likelihood of tracer- 
black hole binary formation, as exemplified by model 5 which 
showed no tracer bound to a black hole at any time during 
the course of the simulation. This may be due to two coun- 
teracting effects; a reduced softening will allow for closer 
encounters, but at the same time will increase the heating 
of the tracer population. This reduces the central density of 
tracers more rapidly than for a system with a larger soft- 
ening and decreases the probability of the close encounters 
which result in the formation of a bound system. 

Given that our simulations were performed using non- 
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Figure 8. Plots showing the decline in tracer number within 
500 pc. Top panel: model 4 (solid), model 5 (short dashed), model 
6 (dotted), model 7 (long dashed). Bottom panel: model 8 (solid), 
model 9 (dotted), model 10 (long dashed) and model 11 (short 
dashed). All models shown in the top panel have black holes of 
mass 10 6 Mq; those depicted in the bottom panel have black holes 
of mass 

10 5.0-5.5 Mq Model 4 has 10 4 tracers and black holes 
of mass 10 6 Mq; model 6 is a repeat of model 4 with 10 5 tracers; 
model 5 is a repeat of model 6 with the softening reduced by a 
factor of 100; model 7 is a repeat of model 6 where the halo scale 
length has been doubled. Model 8 has 10 4 tracers and black holes 
of mass 10 5 5 Mq; model 9 is a repeat of model 8 with 10 5 tracers; 
model 10 is a repeat of model 8 where the halo scale length has 
been doubled; in model 11 the black holes have mass 10 5 Mq. 



zero force softening, it follows that the rate of binary forma- 
tion and the dynamical significance of those binaries that do 
form cannot be properly determined from these models. In 
order to address these issues, we ha ve re-simulated models 
4 and 11 using the NBODY4 code iAarsethl|l999f) on the 
GRA PE-6 special-purpose computer board jMakino et alJ 
Il997l) at the Institute of Astronomy, Cambridge. The results 
are discussed in Section [4.31 below. 

We note that in our simulations using the NBODY2 
code we treat the stars as tracer particles within the po- 
tential of the black hole cluster, thereby allowing improved 
efficiency in integration time. It then follows that the black 
holes do not 'see' the stars; they simply feel the potential 
due to the other black holes that are present in the clus- 
ter. The stars are tracers within this dark matter potential 
and their dynamics are unaffected by the stellar population, 
hence we do not expect to observe tracer-tracer binaries. 
Note also that all particles are point masses with no inter- 
nal structure, and hence the formation of a bound two-body 
system requires a three-body interaction, with the third par- 
ticle carrying away the excess momentum. 



4 DISCUSSION 

4.1 Halo Evolution 

In this section, we use simple dynamical models to explain 
the evolutionary trends of both the black holes and tracers 
in our simulations. There is an initial decrease in the num- 
ber of massive particles within 500 pc, where approximately 
75 percent of the tracers reside at the start of the simula- 
tions. As Fig.^Jshows, this decline occurs on the relaxation 
time scale t r h, and may be explained by the energy exchange 
between massive particles within the half-mass radius. The 
inward flow of heat leads to an expansion of the inner halo 
and a decrease in the central density. This is a manifestation 
of the inverted temperature profile of the NFW model (see 
Fig-0! as compared with the more usually considered Plum- 
mer model. We note that such an inversion is not unique to 
the NFW model, a nd other profile s exist which also display 
such an effect (e.g. lDehnenll20o 3). The analytic dispersion 
shows that the massive particles at a radius of approximately 
r s are initially hotter than those in both the inner and outer 
regions. Following this evolution towards an isothermal ve- 
locity distribution, the subsequent increase in the number of 
massive particles within the same region, due to the well un- 
derstood negative specific heat of self-gravitating systems, 
may then be described as a slow evolution towards core 
collapse, similar to what is found in the Plummer model 
JCohrJll980l) . As was mentione d earlier, this evoluti on was 
also noted in the simulations of lHavashi et alJ (|2003). 

4.2 Stellar Heating Mechanism 

The other main effect observed in our simulations is illus- 
trated in Fig. |H] namely the rapid decline in tracer numbers 
within 500 pc. As the figure shows, reducing the mass of the 
black holes comprising the halo decreases the rate at which 
this heating occurs, indicating that the mechanism for the 
increase in mean stellar kinetic energy is equipartition with 
the massive particles. 

If a system consists of multiple populations with dif- 
ferent particle masses, then encounters will tend to lead to 
the establishment of equipartition of kinetic energies, with 
the more massive particles losing energy to the less mas- 
sive particles over time, assuming that velocity was initially 
independent of mass. As a consequence, the more massive 
particles will gradually lose kinetic energy and move towards 
the system centre, whereas the less massive particles will 
move to larger orbits. The time scale on which equipartition 
occurs is related to the half-mass relaxation time. In order 
to check that the rate of heating of the tracers is consistent 
with kinetic energy equipartition, we calculate the theoreti- 
cal rate of energy transfer between the two populations given 
the ir mean kineti c energies at a specified time. This is given 
by JSpitzerlll987l. eq. 2-60): 

dE t _ ^ f <3\^ 2 m t n b r (E b - E t ) , . 

At \ttJ m b (v? + 1£)3/2 ' ^ Z > 

where m is the mass of an individual particle, n is the aver- 
age number density of the black holes, E is the mean particle 
kinetic energy given by 

Ei = -rmvi , (13) 
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Figure 9. The mean tracer kinetic energy (in 10" 10 TV-body 
units) as a function of time for models 4 (solid) and 6 (dotted). 
Evolution towards equipartition rapidly results in the tracers ac- 
quiring energy from the black holes. 

and r is defined by 

r = 47rG 2 milnA. (14) 

The subscripts t and b refer to tracer and black hole re- 
spectively. Since the halo is observed to evolve on the relax- 
ation time scale at the half-mass radius rh, we take n to be 
the mean number density within rh. The value of A in the 
Coulomb logarithm is given by A = 0.4A where N is the 
total number of particles of both populations within rh. 

Equation 11121 can be recast in a more transparent fash- 
ion by defining e, the energy per unit mass of the relevant 
particle species, and replacing nbvrib by pb- If we consider 
the limit where Eb^> Et, then we find 

^ = 4 (67T) 1 / 2 g f b G 2 p b m b InA . (15) 

For a fixed mass density in the gravitationally dominant 
species, the rate of approach to equipartition increases in 
proportion to the black hole mass and decreases as the cube 
of the velocity of the tracer particles. 

The evolution of tracer energy is shown in Fig. [§] for 
models 4 and 6, where the rate of energy change is observed 
to be steady for approximately 2 Gyr. For model 4, the dif- 
ference in the mean stellar kinetic energy in TV-body units 
is approximately 1.12 x 10~ 10 between 0.2 and 2 Gyr which 
corresponds to 10 A-body time units. The rate of energy 
transfer according to equation 11211 varies as a function of 
time, but is ~ 1.5xl0 -11 at 0.2 Gyr and x 10 -11 at 2 
Gyr. Based on these rates we expect to observe an energy 
transfer in the range (1.1 ~ 1.5) x 10" 10 units, which agrees 
well with the energy difference between these two times. 
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Figure 10. One dimensional (line-of-sight) velocity dispersions of 
model 4 at times T = 0, 10 and 20 Gyr (shown with open circles, 
open triangles and open squares respec tively) overplotted with 
Draco velocity dispersion profile from Wilk inson et al.l J20041) 
(filled circles, with lc errorbars). 



Equation 1151 shows that the rate at which the stars 
gain energy from a population of black holes of mass mb 
and mean spatial density p^ is only weakly dependent on 
the details of the mass distribution. In order to confirm this, 
we have performed a number of simulations (not presented 
in Table in which the halo is represented by a Plum- 
mer sphere whose mean density within 3 stellar core radii 
matches that of the NFW profile in model 4. These simu- 
lations confirm that the stellar heating rates for both cored 
and cusped black hole haloes are almost identical for haloes 
with similar mean density in the region initially occupied by 
the stellar distribution. 

The formulae given above are based on the assumption 
that the change in energy of a test particle is due to the 
cumulative effect of many weak encounters rather than a 
single close encounter. It is therefore necessary to be aware of 
the frequency with which tracers experience close encounters 
with the black holes that may render the use of equation 1121 
invalid. For models 4 and 6, we find that such encounters are 
relatively rare, and the results of our energy calculations 
above confirm that the gradual transfer of energy to the 
tracers by the black holes initially within rh is the main 
mechanism by which energy is transferred between the two 
populations. This may be partly due to the value of the 
softening parameter used in those models which tends to 
suppress large-angle scattering events. 
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Figure 11. Overplot of the tracer density in projection for model 
4 and the Plummcr profile for Draco. Projected densities at times 
T = (open circles), 10 (filled triangles) and 20 (open squares) 
Gyr are shown. The fit to the light distribution of Draco, scaled 
such that the central density matches the simulated density at 20 
Gyr, is given by the dotted line. 

4.3 The Effect of Force Softening 

Hard black hole-black hole binaries may have an important 
impact on the dynamics of the entire halo. Softening in turn 
is a key factor in binary dynamics, from the rate of forma- 
tion to their rate of destruction. We have re-simulated mod- 
els 4 and 11 using the NBODY4 code on a GRAPE-6, with 
the aim of quantifying the effects of a smoothed force-law. 
Each of the five re-simulations were performed using a dif- 
ferent random seed to compute the initial conditions; apart 
from the exclusion of tracers, all other parameters were kept 
the same as the original models. Hard binaries (i.e. those 
with -Chiiil > ma 2 where m is the mean mass and a is 
the velocity dispersion) are identified in these simulations 
and treated using regularisation, to enable a n accurate and 
efficient treatment of their int ernal motion <lAarsethlll999t 
Mik kola fc Aarsethlll99l 199&|). The results confirmed that 
the presence of binaries is a very important factor in deter- 
mining the dynamics of the remainder of the system, with a 
single binary in the central region being sufficient to cause 
significant fluctuations in the central density. 

In the simulations with 10 6 Mq black holes, we observe 
a relatively steady rate of binary formation and destruction. 
These binaries tend to form within the inner 1 kpc region 
and often undergo several companion swaps before being 
destroyed or expelled. Although this has interesting conse- 
quences for the evolution of the inner region of the halo in 
particular, the fact that we do not observe any binaries in 
simulations with 10 5 Mq black holes, coupled with the dis- 




Figure 12. Dependence on the individual black hole mass of the 
time taken for the initial tracer count within 500 pc to drop by a 
factor of 2. The points show this time for models 2-4, 8 and 11; 
the curve shows the fit described by equation 1161 . 



cussion below that black hole masses of greater than 10 5 Mq 
are inconsistent with the preservation of the stellar popula- 
tion in Draco, implies that our overall conclusions based on 
the simulations with softening are unaffected by the inclu- 
sion of softening. 



4.4 Comparison with the Draco dSph 

A comparison of the observed and simulated data can be 
made by means of the velocity dispersion and projected den- 
sity shown in Figs. llUl and llll respect ivelv. The former shows 
the velocity dispersion profile of Draco overplotted on the 
dispersions for model 4 at three sampling times, while the 
latter gives the light distribution at the same times. Also 
plotted in Fig. II II is a Plummer profile, which is known to 
provide a good fit to the light distribution in the inner region 
of Draco. The profile has been scaled such that the central 
density matches that of the core density in the simulation 
at 20 Gyr. The final model density profile in projection is 
clearly too extended to be compatible with Draco's observed 
stellar distribution. 

If the halo of Draco is indeed composed of massive black 
holes, the current extent of the stellar distribution can be 
used to place an upper limit on the mass of individual black 
holes, if we assume that there is no spread in the mass of 
black holes in a halo. In Fig. 1121 we plot the time taken for 
the initial tracer count within 500 pc to drop by a factor of 2 
for models 2-4, 8 and 11. We would expect that this time t 
would be related to £ r h for the models and should therefore 
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scale with mass approximately as M 0,8 , in the mass range 
we are investigating. In fact, the relation 



0.5 



t = 



9.549 x 10 7 



(16) 



with t in Myr and M in M Q , provides a good fit to the re- 
sults from the simulations as shown in Fig. 1121 If we require 
that the process described above takes at least 5 Gyr, we are 
able to set an upper limit on M of 9.7 x 10 4 M©; for 10 Gyr, 
an upper limit of 4.4 x 10 4 Mq would be implied. It should 
be noted that a drop in the central tracer count by a fac- 
tor of 2 often corresponds to an increase in the scale radius 
of the tracer population by a factor of a few and therefore 
these upper limits are probably overestimates. Thus, black 
hole masses of less than 10 5 Mq are required to prevent sig- 
nificant evolution of the tracer density in a Hubble time; 
if Draco's initial stellar distribution was similar to that ob- 
served today then black hole masses in excess of 10 5 M can 
immediately be ruled out. 

It is worth noting that black hole masses below 10 5 M 
are a priori less interesting as dark matter candidates for a 
number of reasons. First, as L085 have shown, black hole 
masses of approximately 10 6 M© are required to explain the 
heating of the Milky Way disk. Second, haloes comprising 
lower mass black holes evolve significantly more slowly than 
haloes of more massive holes and a halo which formed with 
a dark matter cusp would, therefore, not develop a core in 
less than a Hubble time. We therefore conclude that if it 
could be demonstrated that Draco's stellar distribution has 
not changed significantly over a Hubble time, then the two 
primary motivations for adopting black holes as dark matter 
candidates would disappear. 

Another consideration is whether our restricted choice 
of initial conditions would affect our black hole mass limits. 
We have made a number of simplifying assumptions in our 
choice of initial conditions. First we assume spherical sym- 
metry for our halo distributions. This is a reasonable start- 
ing point since dSphs such as Draco are observed to have an 
ellipt icity of the stellar distribution of only e p ~ 0.3 iMated 
1998). The ellipticity of the unde rlying potential is typically 
about 1/3 that of the density fe.g. lBinnev fc Tremainell987l) 
and the spherical symmetry is probably a good approxima- 
tion to Draco. Second, we assume that the dSph evolves in 
isolation - as discussed earlier, observations of Draco sug- 
gest that it has been only weakly affected by the tidal field 
of the Milky Way. Third, we have assumed that all our black 
holes have the same mass - this effectively means that we 
are assuming a strongly peaked black hole mass function. 
The key point with regard to our models is that includ- 
ing any of the above effects e.g. a spectrum of black hole 
masses or an external tidal field w ould serve to speed u p 
the evolution of the black hole halo iGiersz fc Heggidll997t) . 
Our simulations therefore represent a lower limit to the ex- 
pected rate of halo evolution. The inclusion of a spectrum 
of black hol e masses could reduce the time to core collapse 
to ~ Ut llGiersz fc Heggidll996l) . From Table 1, allowing 
for the initial evolution of the NFW profiles towards a core, 
models in which the peak of the black hole mass function 
is less than about 10 5,5 M© are still unlikely to reach core 
collapse within a Hubble time. On the other hand, models 
with a mass peak around 1O 6 M0 or above will almost cer- 
tainly achieve core collapse. This adds further weight to our 
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Figure 13. Top panel: growth of the scale radius of the tracer 
distribution for model 11 as a function of time (solid line). The 
dashed line shows the straight line fit used to describe this growth 
(see text for details). Bottom panel: tracer scale radius evolution 
for models 2 (long-dashed curve), 4 (dotted curve) and 8 (short 
dashed) together with predicted growth rates (solid lines) based 
on a rescaling of the fitted curve for model 11. 



conclusion that black holes of 10 M© and above are unlikely 
dark matter candidates. 



5 COMPACT INITIAL CONDITIONS 

The results of the previous section demonstrate that the 
present stellar distribution of Draco can be preserved for a 
Hubble time only for black hole masses less than 10 5 M©. In 
this section we consider an alternative possibility, namely 
that the stellar distribution was initially significantly more 
compact than it is at present and was subsequently heated 
to its current condition by the processes described above. We 
now present a simple scaling argument which demonstrates 
that the initial conditions required are not unreasonable for 
black hole masses of about 10 5 M Q . 

In order to quantify the evolution of the tracer pop- 
ulatio n, we fit t hree dimensional density profiles of the 
form dZhadll996T) 



p(r) = 



Pari 



(17) 



These profiles incorporate the Plummer profile (a = 0, /3 = 
2, 7 = 5) but they also include outer profiles which fall off 
either faster or slower than r~ 5 . The curves in Fig. 1131 show 
the evolution of the tracer scale radii r c obtained from fits 
of this form to four of our models. 

The top panel of Fig. ll3l shows the evolution of the scale 
radius of the tracer distribution in model 11 during the first 
2 Gyr. Until this time, more than 70 percent of the tracers 
lie within one scale radius r s of the halo. This is important, 
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as we expect the simple scaling arguments given below to 
hold only during the time when the tracers are confined to 
the central regions of the halo, where the density profile is a 
single power-law. The solid line in Fig. If 3 1 shows the straight 
line fit to the temporal increase of r c 

r c (t) =0.28(l + 9.467 (^-)) . (18) 

where r c is in kpc. In the regime where the majority of trac- 
ers lie in the volume where the halo potential goes as 1/r, 
the parameters of the fit should be independent of the mass 
of the black holes and the evolution of the scale radius can 
therefore be estimated from the above equation using the 
appropriate value for t r h- For the other models, this approx- 
imation is valid for less than 2 Gyr as the tracer distribution 
rapidly expands into the outer regions of the halo. However, 
as the bottom panel of Fig. 1131 demonstrates, this simple 
scaling argument yields a reasonable estimate of the initial 
growth rate of the tracer scale radius. At later times, the 
heating rate falls below the estimated value as the tracers 
no longer lie solely in the single power-law region of the halo. 

We can extrapolate the scale radius evolution estimated 
from equation 1181 to consider the evolution of tracer popu- 
lations whose initial distribution is more compact than those 
considered so far. This extrapolation backwards in time is 
valid unless the initial r c is so small that the stellar popula- 
tion would become self-gravitating. For example, for model 
2, if we require that after 10 Gyr we have r c — 232 pc, then 
this implies an initial radius of 1 pc for the stellar popula- 
tion, which corresponds to a dense and self-gravitating stel- 
lar cluster. However, for model 11 the required initial ra- 
dius is 48 pc. If this population resided in the centre of the 
halo, then the halo mass within its volume would be com- 
parable to the stellar mass which would not, therefore, be 
self-gravitating. This suggests that the heating of the stel- 
lar population by a halo of 10 5 Mq black holes need not be 
catastrophic for a dSph such as Draco, provided that its ini- 
tial stellar population is sufficiently compact. It is less clear 
whether the same holds for black hole masses of 10 6 Mq. In 
this case, the initial cluster has scale radius 8pc and since it 
would therefore be almost self-gravitating, the simple scaling 
argument presented here breaks down. 

It is important to emphasise at this point that the heat- 
ing rate depends on the relative velocities of the black holes 
and tracers. A more compact, but not self-gravitating, initial 
tracer density distribution has a lower velocity dispersion 
and hence the rate at which it is heated is initially greater 
than we have estimated here. However, this increases the 
initial heating rate by at most a factor of two for the range 
of initial conditions we consider and so will not substantially 
affect our conclusions. There are additional heating effects to 
be considered for an initially small system, such as relaxation 
due to vN fluctuations in the number of black holes within 
the stellar distribution. These would make early expansion 
proceed more rapidly than otherwise. While it might be pos- 
sible to balance this effect to some extent by the stabilising 
effects of self-gravity, we still conclude that constructing a 
dynamical history for Draco which results in a present-day 
stellar distribution with a half-light radius of 232 pc may be 
problematic for the scenario considered here. Further work 
is needed to establish the exact degree of fine tuning of the 
initial conditions which would be required. 



If the stellar population of Draco began in a compact 
state and its current extent is the result of heating by black 
holes, it is natural to ask whether this process would lead 
to any clear observational signatu res in either th e stellar lu- 
minosity or velocity distributions. ISpitzerl il987f) has shown 
that weak encounters in iV-body systems produce a pop- 
ulation of weakly bound objects with a power-law density 
distribution. The logarithmic slope of this density profile de- 
pends on the gravitational potential, which is assumed also 
to be a single power-law. L085 showed that strong encoun- 
ters generate a profile with a different logarithmic slope to 
that generated by weak encounters. In our simulations, how- 
ever, both weak and strong encounters with the black holes 
contribute to the heating of the stellar distribution. As a re- 
sult, there is no simple expression for the outer slope of the 
resultant light distribution. If the tail generated by close 
encounters were sufficiently populated by stars, and could 
therefore be identified in the light distribution of Draco, it 
would argue strongly that the stellar population had been 
heated through close encounters. On the other hand, our re- 
sults demonstrate that the absence of such a tail does not 
preclude the black hole model. For example, in model 5 the 
outer profile of the stellar distribution after 10 Gyr goes ap- 
proximately as r -4,6 which would be indistinguishable from 
a Plummer law, given the magnitude of t he error bars on the 
observ ed luminosity at such radii (e.g. lOdenkirchen et ail 
l2001al) . We note that the change in the slope of the light 
distribution in the outer parts of Draco r ecently identified 
from deep imaging JWilkinson et al.ll2004h might be indica- 
tive of such a power law halo. 

The velocity distribution which arises in the black hole 
halo scenario is also difficult to identify unambiguously from 
observations. One would expect to observe a significant pop- 
ulation of stars on radial orbits in the outer regions of a 
heated model as a result of the heating process. Indeed, af- 
ter 10 Gyr, the anisotropy of the velocity distribution in 
model 5 is strongly radial in the outer parts. However, the 
corresponding projected velocity distribution in that model 
does not display any features which would unambiguously 
point to the strong radial bias in the velocity distribution. 
We conclude that the simulations presented here do not pre- 
dict a unique observable signature of heating by black holes. 
In order to address this rather unsatisfactory situation, wc 
will investigate the evolution of initially compact dSphs in 
a future paper. 

In our simulations we have shown that a halo which 
is composed of massive black holes and which initially has 
an NFW density profile can be transformed into a cored 
halo through dynamical evolution. A number of authors have 
considered the role of small numbers of black holes in trans- 
forming th e density profiles of cusped stellar systems into 
cores (e.g. lHemsendorl| | 2003|) or cusped haloes into cored 
haloes (e.g. iMerritt &; Milosavlievial2002r) . In these papers, 
the black holes were a minor constituent of the mass distri- 
bution, contrary to the case we consider in which the black 
holes are the dark matter and constitute the bulk of the 
gravitating mass. We are not aware of any simulation of 
the evolution of a compact stellar system in the presence of 
an extended halo of black h oles. Of more direct re l evance 
to our results, we note that iHernandez fc Gilmord (Il998ft 
concluded that black holes with masses above 10 s Mq could 
not constitute more than 1/8 of the total mass of dwarf disk 
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galaxy haloes, as in this case dynamical friction would lead 
to a concentration of mass at the centre of the dwarf and 
a consequent observable signature in the rotation curve of 
the dwa rf. Our results cannot be compared directly with 
those of iHernandez fc Gilmord as their calculation assumes 
a smooth background halo distribution with a uniform den- 
sity core, while in our models no su ch smooth background 
exists. However, it is interesting that IHernandez fc Gilmorel 
came to similar conclusions about the upper limit to the 
mass of any putative black hole population. 

We conclude this section with a discussion of potential 
criticisms of our black hole halo model. First, one might 
argue that the probability that we would observe a signifi- 
cant number of dSphs at a similar phase of their evolution 
should be very small, which would render our compact initial 
conditions scenario implausible. This objection can be coun- 
tered if we assume that dwarf elliptical galaxies and dSphs 
were originally members of the same class. Dwarf ellipticals 
have high stellar densities compared with the low densities 
of the dSphs, and observations have produced evidence for 
the presence of black holes in these galaxies with masses of a 
few million solar masses in a number of these systems (e.g. 
M32: IVerolme et all 1200^ . As we have seen, heating of a 
compact stellar system by a black hole halo could produce 
a low surface brightness object similar to a dSph. On the 
other hand, in some cases a single black hole might become 
embedded in the compact cluster protecting it from further 
shredding and causing the stellar density to remain high, as 
in M32. If there is a characteristic black hole mass, then the 
evolutionary time scales of all the dSphs would be similar 
and hence it is not surprising that we observe them today 
in similar evolutionary states. 

A second objection to the hypothesis that black holes 
constitute the dark matter in dSph galaxy haloes arises from 
observations of the Fornax dSph. Fornax, the most luminous 
Local Group dSph, is surrounded by a population of five 
old, globular clusters. The presence of these clusters in an 
extended distribution about the galaxy is difficult to recon- 
cile with the presence of a large number of black holes whose 
masses are comparable to those of the clusters, because close 
encounters between the clusters and the black holes would 
be expected to lead to the rapid disruption of the clusters. 
While this is indeed perplexing, we note that if the dark 
matter is composed of much lower mass particles (e.g. sub- 
atomic particles), then the extended radial distribution of 
the Fornax clusters is equally surprising, as dynamical fric- 
tion should draw the clusters into th e centre of the galaxy 
in a fraction o f the Hubble time (e.g. iTremaine et alJll975t 
lOh et all2000h . Thus the origin of the Fornax globular clus- 
ters remains an open question and, in the absence of further 
data, does not provide a means of distinguishing between 
the conventional model for the dark matter and the pro- 
posal considered here. 

There are also two sets of objects in the Milky Way halo, 
observations of which appear to argue against black holes as 
major constituents of the dark matter. If black holes were 
ruled out for the Milky Way halo, it would make it less likely 
that the dark matter in dSph galaxies could be composed 
of black holes, since disruption of dSphs over the lifetime of 
the Milky Way would pollute its halo with black holes. First, 
deep imaging of the Milky Way globula r cluster Pal 5 has 
revealed the presence of long tidal tails llOdenkirchen et alJ 



I2001bl.l2003h extending to distances of about 2 kpc from the 
cluster. Estimates of the drift rate of stars along these tails 
suggest that the time required for the tails to reach their 
obser ved extent is approximately 2 Gyr llOdenkirchen et alJ 
2003). It has been pointed out that the spatial integrity 
of the streams, which implies that they are kinematically 
cold, may place constraints on the properties of any massive 
substructures prese nt in the halo of the Milky Way (e.g. 
iDehnen et all I2004T) . In particular, simulations by Moore 
(priv. comm.) suggest that the tails of Pal 5 rule out a signif- 
icant population of black holes with masses of about 10 6 M 
in the halo of the Milky Way. However, we note that simi- 
lar arguments apply to the substructure which is ubiquitous 
in standard cold dark matter simulations and thus the tails 
of Pal 5 do not constitute a clear obs ervational test of the 
black hole model. We also note that IDehnen et alJ J2004T) 
have shown that iV-body simulations of the tidal disruption 
of Pal 5 in a smooth Milky Way halo are unable to repro- 
duce the observed internal structure in the stellar density 
distribution along the tails. Those authors suggest that the 
presence of an amount of halo substructure may be required 
to account for the non-uniform distribution of stars along 
the tails of Pal 5. 

Recent work on the distribution of separations of wide 
stellar binaries in the halo of the Milky Way may con- 
stitute the strongest observational argument against the 
black hole halo mode l. In a comprehensive study of this is- 
sue, lYoo et all i2004f) show that encounters between wide 
binaries and massive compact objects would result in the 
power-law tail of binary separations becoming steeper with 
time. They argue that the absence of a cut-off in the distri- 
bution of wide binary separations in a large sample of nearby 
halo binaries rules out the presence of compact objects with 
masses above 43 Mq at the standard local halo density. This 
analysis assumes that all the binaries in their sample have 
been part of the Galactic halo for at least 10 Gyr and during 
that time have experienced a density of dark matter compa- 
rable to the local halo density. Given that most of the halo 
binaries were probably formed in star clusters, it is not clear 
how long they have been exposed directly to potentially dis- 
rupting encounters with halo substructure. Further, if the 
binaries originated in halo star clusters on highly elongated 
orbits, they may have spent much of their lives in regions 
where the dark matter density is significantly lower than in 
the solar neighbourhood. We conclude that without infor- 
mation about the actual distribution of ages for the halo 
binaries (i.e. how long they have been part of the halo) and 
knowledge of the orbital distribution of the binaries in the 
Galactic halo, the population of wide binaries is not yet suf- 
ficient to rule out definitively massive objects as a significant 
component of the dark matter. However, these data have the 
potential to rule out black hole dark matter and therefore 
merit further analysis to determine their orbital properties. 



6 SUMMARY AND CONCLUSIONS 

The high mass-to-light ratios of certain Local Group dwarf 
galaxies suggest that they are among the most dark matter- 
dominated stellar systems known in the Universe. This 
makes them ideal laboratories in which to investigate the 
nature and properties of dark matter. In addition, their rel- 
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ative proximity has allowed the accumulation of a wealth 
of high quality stellar kinematic data in recent years. We 
present here the results of simulations where the modelled 
system is assumed to consist of a dwarf galaxy-sized stel- 
lar population, which initially resides in the centre of an 
extended, massive halo with an NFW profile. The halo par- 
ticles are assigned masses of between 10 5 Mq and 10 7 Mq 
and are modelled as massive black holes. The mass of the 
system is chosen to match that of the Draco dSph, a Galactic 
satellite for which evidence has recently bee n presented of 
the e xistence of a massive dark matter halo llKlevna et alJ 
l200lT) . The stellar distr ibution appears to have suffered only 
minor tidal distortion jKlessen et all2003t fwilkinson etlll 
12004) . This, together with the availabilit y of stellar velocity 
measurements up to four Plummer radii na et al 12002: 
IWilkinson et al]l2004h . makes Draco a useful test candidate 
for comparison with simulations. 

The first question which our simulations address is 
whether the black hole population evolves too quickly to 
be consistent with the potential observed for Draco. Perusal 
of the very short central relaxation times (Table0 indicates 
that we might have expected the systems to undergo core 
collapse and re-expansion within a Hubble time, with sub- 
sequent dissolution in a Galactic tidal field, had the initial 
cond itions been closer to a Plummer model llLee fc Ostrikerl 
Il987f) . However, due to the inverted initial temperature pro- 
file of the NFW model, the evolution in our simulations 
is slower and initially in the opposite direction, towards a 
lower central density with a cored rather than a cusped den- 
sity profile. Only systems with black hole mass greater than 
10 6 ' 5 Mq appear to evolve too rapidly towards core collapse. 
For the lower end of the mass range investigated, evolu- 
tion leads to a profile with a density profile that is more 
consistent with observations of dSphs jKlev na et al1 l2003: 
lMagorriarJl2003h than the initial NFW model. 

A more serious set of issues is raised by the evolution of 
the tracer populations. We find that in all of our dSph-sized 
models, the central tracer density shows a very rapid decline 
in the first Gyr, corresponding to over-efficient heating by 
the black holes. This effect can also be observed via the 
increase in the tracer velocity dispersion, which at 10 Gyr 
cannot be reconciled with the observed Draco dispersion. 
If the initial stellar distribution is chosen to match that of 
Draco, then only for black holes of mass below about 10 5 Mq 
is the heating rate reduced to a level which could perhaps 
be compatible with the observations. We conclude that the 
haloes of dSph galaxies cannot be composed of black holes 
more massive than 10 5 Mq if their initial stellar distributions 
were similar to those observed at the present time. 

We have also performed simulations in which the scale 
length was set to twice that used to model the halo of Draco. 
The heating of the stellar population in these cases was 
somewhat reduced, and more than 10 percent of the stel- 
lar population which initially occupied the central 500 pc of 
the halo remained in this region after 10 Gyr for black hole 
masses of 10 6 Mq , as opposed to a loss of greater than 95 
percent for the models with smaller r s . However, the veloc- 
ity dispersion of the evolved tracers within the same region 
is still uncomfortably large. Notwithstanding the question 
of how to suppress the heating of the stellar population, the 
evolution of the haloes of these models is potentially very 
interesting. As Fig. [3]shows, the inner regions of the haloes 



evolve towards a more cored profile. Thus, these models may 
be able to explain the observations of L SB galaxies which 
appe ar to have large dark matter cores <lde Blok fc Bosmal 
120021) . 

By constraining the amount of evolution allowed to take 
place for the stellar population within a given time scale, we 
are able to place a new upper limit on the mass of indi- 
vidual black holes which constitute the dark matter halo 
of Draco. If we require that the stellar count within 500 pc 
of the halo origin decreases by no more than a factor of 2 
within 5 Gyr, we obtain an upper mass constraint of ap- 
proximately 10 s Mq; the equivalent condition required over 
a time scale of 10 Gyr sets a limit of 10 4 ' 6 Mq . These num- 
bers call into question the hypothesis that massive black 
holes could be a significant component of the dark matter: 
black holes with masses below 10 5 Mq are less interesting as 
dark matter candidates because they cannot explain either 
the heating of the Milky Way disk (L085) or the origin of 
cores in dark haloes (Section 14.41 . Therefore, we have also 
investigated the possibility that a dSph whose initial stellar 
distribution was more compact than those of the present- 
day dSphs could evolve into an object resembling Draco as 
a result of heating. We have presented a simple scaling ar- 
gument which implies that if the black holes have masses of 
10 5 Mq then an initial stellar distribution with scale radius 
48 pc would evolve into an object of similar extent to that 
of Draco in 10 Gyr; such a scenario may permit models with 
black hole masses somewhat larger than 10 5 Mq and cor- 
respondingly more compact initial conditions. However, we 
note that this scenario may require very finely tuned initial 
conditions because y~N noise in the black hole distribution 
may inflate the heating rate at early times. 

In summary, we find that scenarios in which the dark 
matter is in the form of massive black holes, with masses 
between 10 5 Mq and a few xlO 5 Mq, may provide a viable 
channel for the production of cored, rather than cusped, 
haloes but only if the initial stellar distributions of the dSphs 
were considerably more compact than the currently observed 
ones. Further simulations will be required to determine if 
this picture is, in all details, consistent with the observed 
data on the Local Group dSphs. 
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